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Thermodynamic perturbation theory is employed to derive analytical expressions for the equi- 
librium linear susceptibility and specific heat of lattices of anisotropic classical spins weakly coupled 
by the dipole-dipole interaction. The calculation is carried out to the second order in the coupling 
constant over the temperature, while the single-spin anisotropy is treated exactly. The temperature 
range of applicability of the results is, for weak anisotropy (A/ksT -C 1), similar to that of ordinary 
high-temperature expansions, but for moderately and strongly anisotropic spins (A/ksT > 1) it can 
extend down to the temperatures where the superparamagnetic blocking takes place (A/ksT ~ 25), 
provided only the interaction strength is weak enough. Besides, taking exactly the anisotropy into 
account, the results describe as particular cases the effects of the interactions on isotropic (A = 0) 
as well as strongly anisotropic (\A\ — > 00) systems (discrete orientation model and plane rotators). 
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I. INTRODUCTION 



In paramagnetic systems the relative weakness of the dipole-dipole interaction between magnetic ions results in 
characteristic temperatures lying in the range of 0.01-0.1 K. Besides, the dipolar coupling usually coexists with the 
i-^j \ exchange interaction. However, for superparamagnets (nanoscale solids or clusters whose net spin rotates thermally 
<~j . activated in the magnetic anisotropy potential) , the exchange or other competing interactions can usually be discarded 
and the dipole-dipole interaction can be studied in pure form. In addition, the size of their typical magnetic moments 
(S ~ 10 2 -10 5 ) shifts the relevant temperatures up to the range of a few Kelvin, facilitating greatly the experimental 
study of this interaction.El 

The calculation of the relevant statistical-mechanical quantities constitutes a formidable problem in most many-body 
^ ■ systems. Apart from various specific solution ansatzs, a number of systematic expansions (in the density, coupling 
parameter.,, .etc.) have been developed for weak interactions. a In spin and dipole systems, the moment method of 
Van Vleckd'Q (a high-temperature expansion of the partition function) permits to study the equilibrium properties in 
the absence of cooperative phenomena. This technique is one of the few analytical tools available to handle systems 
coupled by the dipole-dipole interaction, due to the long-range and reduced symmetry of this interaction. 
' An important property of superparamagnets is their magnetic anisotropy, which.. results in a number of spin ori- 
entations of minimum energy separated by potential barriers. For uniaxial spinsO the characteristic time for the 
■ thermo-activated rotation of the spin over the anisotropy barrier A can approximately be written as 
Ctf ■ 

g: r~T exp(A/fc B r), (1.1) 

1 ■ 

where To is weakly temperature dependent and takes values To ~ 10~ 10 -10~ 12 s for magnetic nanoparticles. Then, for 
a given measurement time, t m , the system exhibits its thermal-equilibrium response when the condition of superpara- 
magnetism, t m 3> t, is obeyed, which corresponds to the temperature range (in units of A/ks)'. 



Ht m /r Q ) > A/k B T > 0. (1.2) 



In "static" measurements (i m ~ 1-100 s), due to the small value of To, this equilibrium range extends down to very 
low temperatures (25 > A/HbT), showing that the naive ascription of superparamagnetism to the range in which 
"the thermal energy is comparable or larger than the anisotropy energy" (1 > A/k^,T) is unduly restrictive. Indeed, 
as T decreases the system displays different behaviors ranging from almost isotropic (A/k^T <§; 1), then moderately 
anisotropic {A/ksT ~ 1), and eventually strongly anisotropic {A/ksT 3> 1) without leaving the equilibrium regime. 
Therefore, descriptions based on the assumption of isotropic behavior or the opposite discrete-orientation or plane- 
rotator approximations (for easy-axis and easy-plane anisotropy), necessarily have a reduced range of validity in 
superparamagnets. 

Due to the mentioned characteristics of the dipole-dipole coupling and the difficulties introduced by the anisotropy, 
most rigorous calculations in interacting superparamagnets have been done by numerical simulation techniques In 
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this article we apply thermodynamic perturbation theoryu to calculate analytically the linear susceptibility and the 
specific heat of lattices of uniaxial classical spins coupled by the dipole-dipole interaction, accounting exactly (non- 
perturbatively) for their anisotropy energy. Along with the study of the. dependence on the shape of the system of 
certain quantities (due to the long-range of the dipole-dipole interactions), our treatment permits to investigate the 
effects of the strength and sign of the anisotropy, as well as of the orientational distribution of anisotropy axes. 

We find that for systems with axes oriented at random the corrections to the specific heat and the linear suscepti- 
bility become independent of the anisotropy (at least to second order in the interaction coupling) in certain spatial 
arrangements of the spins (e.g., cubic or completely disordered). The latter is a generalization to interacting systems of 
the well known absence of anisotropy effects on the equilibrium linear susceptibility of dipoles with random anisotropy 
(discussed in Ref. |i|). However, apart from the important exception of random axes, the anisotropy is an essential 
element in the determination of the corrections due to the interactions. This is illustrated with the response of systems 
with parallel anisotropy axes, where an ordinary high-temperature expansion, either disregarding the anisotropy or 
in the discrete orientation limit, poorly describes the susceptibility curves (computed for comparison by Monte Carlo 
simulation), while the thermodynamic perturbation theory describes the results with reasonable accuracy. 



II. THERMODYNAMIC PERTURBATION THEORY FOR INTERACTING DIPOLES 



In this section we introduce the spin system studied and discuss the application of perturbation theory to calculate 
approximately thermodynamic quantities. 



A. Hamiltonian of a system of interacting anisotropic spins 

Let us consider a system of N magneto-anisotropic spins coupled by the dipole-dipole interaction. The magnetic 
anisotropy energy is assumed to be uniaxial 

E a = -Aj2(si-ni) 2 , (2.1) 

i 

where A is the anisotropy parameter (for magnetic nanoparticles A = KV, where K and V are the anisotropy constant 
and volume), and Sj and fU, are, respectively, unit vectors along the magnetic moment and anisotropy axis of the ith 
spin. The dipole-dipole interaction energy can be written as 

i>j 

where m is the magnitude of the magnetic moment, a is an appropriate characteristic length (see below), and 

G « = 3" (Zrijfij - 1), (2.3) 
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r J' ?ij ~ r ij/ r ij ■ (2-4) 



Here 1 is the unit tensor and 7% the vector joining the sites i and j (measured in units of a). The action of a tensor 
dyadic T = uv on a vector w is the usual one (uv)w = u(v • w), and hence the tensor G^, when multiplied with sj, 
gives (except for a constant) the field at the position of the zth dipole created by Sj. 

For notational simplicity we are assuming that the parameters characterizing the different dipoles are identical, 
but it is immediate to generalize the expressions for different anisotropy constants, magnetic moments, volumes, etc. 
Similarly, although the assumption of uniaxial anisotropy is not necessary until the end of the calculation we made it 
here for definiteness. Concerning the characteristic length a, it is defined in such a way that a 3 is the mean volume 
around each spin. Thus, a is the lattice constant in a simple cubic arrangement, and for nanoparticles the volume 
concentration of particles is c = V/o? . 

Finally, introducing (3 = l/k^T and the following dimcnsionlcss quantities (anisotropy and coupling constant 
relative to the thermal energy) 

cr= A ^ = fi m 2 1 

k^T ' 47ra 3 k^T ' 



2 



we can write the total energy E = E a + as 

- (3E = aJ2& ■ Atf + ZdY, 10 ^ ■ ( 2 - 6 ) 

i i>j 

Note that the interaction strength can also be measured by the temperature independent coupling parameter 

h d = £ d /2a , (2.7) 

which is the magnitude of the field (measured in units of the maximum anisotropy field uqHk = 2A/m) produced at 
a given position by a dipole at a distance an 



B. Equilibrium linear susceptibility and specific heat 

The thermal-equilibrium average of any quantity B(s±, . . . , sn) is given by 

{B) = \ J dTB ex P("^) ' ( 2 - 8 ) 

where Z = J dT exp(— j3E) is the partition function. In classical spins the different states correspond to different spin 
orientations so that dT = Y[ { dQ,i, with dfli — d 2 Si/2ir. 

The linear susceptibility is defined as the derivative of the magnetization jr(Y^i^i) with respect to the external 
magnetic field (which is the experimentally manipulable quantity in contrast to the internal macroscopic field). 
However, from basic statistical mechanics we know that the response to a probing field AH can be obtained in terms 
of suitable averages of the net spin taken in the absence of AH. If in addition there is no external bias field applied, 
the susceptibility is simply given by 



Horn 2 1 



«, = V(Si-/i), (2.9) 



k B T N x zn ^ 

where h is a unit vector along AH and s z is the field projection of the net moment. The specific heat at constant 
volume c v = d(E)/dT can be obtained directly from Z as 

g-^O^-^O.*). (2.10) 

where to take the cr-derivative, the coupling parameter £d is expressed as £d = 2ahd [Eq. (|2.7| )1. As in the calculation 
of X: we only consider the zero-field specific heat.113 



C. Thermodynamic perturbation theory 

We shall now use thermodynamic perturbation theory}] to expand the Boltzmann distribution W = Z^ 1 exp(—f3E) 
in powers of £d ■ This will lead to an expression of the form 

W = W a (1 + CdFi + \(lF 2 + ■■■) , (2.11) 

where F\ is linear in E& (and hence quadratic in the spins), F% is up to quadratic in E& (quartic in the s*j), and 

iy a = Z- 1 exp(-/3S a ) , (2.12) 

is the Boltzmann distribution of the noninteracting ensemble. Therefore, the calculation of an observable (B) is 
reduced to the calculation of averages weighted by W a (denoted (-) a ) of typically low grade powers of the spin 



variables: {B) a , (BFi) a , (BF2) a , etc. An ordinary high-temperature expansion corresponds to expand Eq. ( 2.12| ) 
further in powers of j3 = l/k-oT. In the present calculation, however, the averages are kept weighted over W a , so they 
will be exact in the magnetic anisotropy and only perturbational in the dipolar interaction.EU 

A convenient way of performing the expansion in powers of £d is to introduce the Mayer functions fij defined by 
1 + fij = exp(^dWy ), which permits to write the exponential in the Boltzmann factor as 
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cxp(-l3E) = cxp(-f3E a ) + /y). (2.13) 

i>j 

Expanding the product to second order in the gives 

+ fij) = 1 + CdGx + \&G 2 + 0(3), (2.14) 

i>j 

whereQ 

Gi=2w tfl (2.15) 

G 2 = ^ W % + X X] u) i3 u} kiqik:jlQa:jk, (2-16) 

and the symbol gjfe^/ annihilates terms containing duplicate pairs: qik-.j i = 1 (2 — (5,;fc — ^;)(1 + 5jfc)(l + 

To obtain the average of any quantity £? we introduce the expansion (2.14) in both the numerator and denominator 
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of (B) = JdTB exp(— @E)J JdT exp(—f3E), and work out the expansion of the quotient, getting 

(B) — (B) a + £d [ (B Gi) a — (B) a (Gi) a ] 
+ ie d {(BG 2 ) a -(B) a (G 2 ) a 

-2(Gi) a [{BGi) a -(B) a (Gi) a ]}. 

However, since in our case W a (— Si) = W a (si) (because the single-spin anisotropy has inversion symmetry and there 
is no bias field) and in G\ = J2i>j u ij a dipole does not interact with itself, the result (Gi) a = holds. Under these 
conditions we finally find the simpler form 

(B) ~{B) a + t d {BG 1 ) a 

+ ^[(BG 2 ) a -(B) a (G 2 ) a ], (2.17) 

To complete the calculation we need to obtain averages of low grade powers of s weighted by the noninteracting 
distribution (moments), which is the only place where one needs to specify the form of E a . We can write the 
susceptibility and s peci fic hea t to second order in £d using up to fourth order moments, which are calculated in 
Appendix [A] [Eqs. (Al) and (|A2|)] for the uniaxial distribution W a cx exp[cr(s • ft) 2 ]. For instance, the first two 



moments read ((ci • s)) a = and 

1 - S 2 

((ci • s)(c 2 • s)) a = — - — ci • c 2 + S 2 (ci • n)(c 2 ■ n) , 

where the c„ are arbitrary constant vectors and Si is the average of the Ith Legendre polynomial Pi(z) over W a ; 

S l (a) = (P l (s-n)) a . (2.18) 

In particular, So = 1 and S^cr) = \ (3(s • n) 2 — 1^ can be expressed in terms of familiar special functions (error 
functions or the Dawson integral), while oth er Si appearing in higher order moments, can be obtained by means of 
the recurrence relation they satisfy [Eq. (|A8[)]. 



III. ANALYTICAL EXPRESSIONS FOR THE SUSCEPTIBILITY AND SPECIFIC HEAT 

We shall now present the general form of the perturbative expressions for the susceptibility and specific heat and 
discuss their properties in some important cases. 
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A. Linear susceptibility 



Using the results of the previous section, we have averaged the square of the field projection of the net spin 
\B = (y]j Si ■ h) 2 }, which yields an expansion for the equilibrium linear susceptibility [Eq. (|2.9| )1 of the form 

X=^(aa+Ua 1 + i 2 e d a 2 ) . (3.1) 

The lengthy general expressions for the coefficients a n , which include sums over the lattice of the dipolar tensor 
and the Si, are written in full in Appendix |b| 

The coefficients a n simplify notably for some orientational distributions of the anisotropy axes. For systems with 
parallel axes (e.g., single crystals of magnetic molecular clusters, or a ferrofluid frozen in a strong field), the coefficients 
for the longitudinal response read 

00,11 = 3 ^ > 

"1.11 = „ C (3.3) 



4S 2 + 4Sf 



1 iT1^T^ 2 

2 a 2,ll = - 



27 

x [(l-5 2 ) (K-S) + 3S 2 (T-W)] 

7 + 105*2 - 355| + 185 4 
+ 315 

x [(l-S 2 )V + 3S 2 (T-in)] , (3.4) 

where C, 1Z (7£), S, T, U, and V are certain lattice sums whose properties are discussed below. 

To obtain the susceptibility when the aniso tropy axes are distributed at random, we average the general expressions 
for the a n over n, with help from Eqs. ( All ) and ( A12 ) (with s— * n), getting 

flO.ran = ^ (3.5) 
ai,ran = ^C (3.6) 



l«2 iran - ~ (n - S) + 1 (1 - Sl) V . ( 3.7) 



Note that in the limit of isotropic spins (where Si — > 0) the results for coherent axes and for random anisotropy duly 
coincide. 



B. Specific heat 



To obtain the specific heat w e can expa nd directly the partition function in powers of £d by introducing the expanded 
Boltzmann factor [Eqs. ( 2.13 ) and ( 2.14 )] in the definition of Z 



Z = J dTeM-PE a )(l + UGi + ^ 2 d G2) 



z a (i + kj (g 2 > j 



(3.8) 



where the linear term vanishes because (Gi) a = 0. Then c v is obtained [Eq. (2.10)] by differentiating the logarithm 
of the above expansion, which poses no problem of the type of expanding a quotient, since to second order in £d we 
can use ln(l + xt£) ~ xt; 2 . The result has the form 



= o*bo + ^2 



Nk B 



(3.9) 



where the zeroth order coefficient 
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5 " = ^(18S 4 -35S 2 2 + 105 2 + 7) , (3.10) 
315 

gives the specific heat in the absence of interactions. 

The general formula for b 2 is given in Appendix |c| [Eq. (CI)]. Again, it simplifies for coherent axes and for random 



anisotropy. In the first case (rij = ft, Vi) we obtain 

&a,|| = 5f - 4aS 2 S' 2 - o- 2 [S 2 S' 2 ' + {S' 2 f])n 

+i(25 2 (l-S 2 )+4aS 2 (l-25 2 ) 

+ ( x 2 {S 2 '-2[S 2 S 2 ' + (S 2 ) 2 ]})v 

+ {sf + AcrS 2 S' 2 + <t 2 [S 2 S 2 ' + (5 2 ) 2 ]}r , (3.11) 

wher e /' = df /da. For randomly distributed axes, on averaging the general expression for 6 2 over n by means of 
Eq. ( |Al3j ), one gets 

&2,ran = (3.12) 

This is the same correction term as that obtained for isotropic spins by WalleiQ and Van Vleck.tl 

C. The lattice sums 

An essential element of the expressions derived for \ an d c v are the following "lattice sums" 

^=§EE^ (3-14) 

n ( 3 - 15 ) 

5 = ^EEE^ G ^- G ^^ ( 3 - 16 ) 

r = ^EE(^ G ^) 2 ( 3 - 17 ) 

w = jf E E • G y • ^ • G j fc • ( 3 - 18 ) 

V = ^EE^- G ^ ( 3 - 19 ) 

i j^i 

(replace h by n in t he formulas for c v ). Considering the structure of these sums along with the form of the dipolar, 
tensor Gjj [Eq. (2^2)], some physical interpretation can be provided for the different terms in the perturbative series.E3 



The first order term a\ incorporates through C the direct action (j — > i) on each spin of the remainder spins, when 
aligned along the probing field. No term of this type appears in the specific heat since the absence of any external 
field yields b\ = 0. The second order terms a 2 and b 2 involve lattice sums including products of with Gjk so they 
take into account the action on a given spin of the others but through intermediate spins {indirect action k — ► j — ► i). 
In particular, if k = i we have the reaction on the ith spin of its direct action on the remainder spins. 

In the next section we shall compute \ and c v for "sufficiently isotropic" lattices, in the sense of fulfilling J2( r x) n = 
12( r v) n — XX r z)"j e -S-i cubic and completely disordered lattices (incidentally, the type of arrangements for which in 
the classical Lorentz cavity-field calculation the contribution of the dipoles inside the "small sphere" vanishes). In 
these lattices we have two important results: (i) 1Z coincides with the more familiar lattice sum 1Z (which justifies the 
notation) and (ii) V = 0. 
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D. Ordinary high-temperature expansions and other approximations 



The effects of the anisotropy are included exactly in our formulas through the anisotropy-weighted averages of the 
Legendre polynomials S2 (f) and S4 (a) , which modulate the different contributions and introduce an extra dependence 
on the temperature via a — A/k^T '. This reflects the fact that the dipolar field at a given site is different if the source 
spins are almost freely rotating (high T, Si — * 0) or, for example, lay almost parallel to their respective anisotropy 
axes (low T, Si — > 1). _ . 

The approximate behavior of S2 and S4 for weak (\cr\ <C 1) and strong (\a\ I) anisotropy isli3 

■■ H«i 

S 2 (a) = { l-^-4^ + "-^» 1 > (3- 2 °) 

• cr < -f 




\a\ « f 

= <! 1 ~ f + 3^ + ' ' ' o- > 1 • (3-21) 
f + ^) + -" a«-f 

The results of an ordinary high temperature expansion, in which all terms in exp(— (3E) are expanded, correspond 
to replace S2 and S4 in the coefficients a n and b n by their weak anisotropy approximations. Besides, taking the 
(7 — ► limit (where Si = 0) we get the known results for isotropic spin s (in a slightly more general form, since the 
terms including V are usually omitted due to the lattices assumedliMa) . In addition, substituting the above strong 
anisotropy formulas in x and c v we get these quantities in the discrete-orientation and plane-rotator cases (with the 
corresponding corrections in powers of l/cr). 



IV. BEHAVIOR OF THE SUSCEPTIBILITY AND THE SPECIFIC HEAT 



In this section we study the features of the susceptibility and specific heat emerging from the analytical expressions 
derived. We shall discuss the shape dependence of these quantities, investigate the dependence on the anisotropy 
(on both its strength and the axes distribution), and finally estimate the limits of validity of the expansions. For 
concreteness, we shall consider the behavior of spins with easy-axis anisotropy (A > 0) in simple cubic lattices. 



A. Shape dependence 

Here we restrict our attention to systems with ellipsoidal shape, which is the geometry usually consider in studies 
of the dipole-dipole interaction because: (i) in the continuous limit the demagnetizing field is spatially homogeneous 
and parallel to the external field when this is along one of the three principal axescll and (ii) it covers as limit cases 
important geometries such as those of disks or long cylinders. 

The long range of the dipolc-dipole interaction leads to a shape dependence of the physical quantities in an external 
fieldu and hence of the linear susceptibility (which is a field derivative). In the expressions obtained, this shape 
dependence is borne by the slowly convergent lattice sums C, S, and U. For the class of sufficiently isotropic lattices 
mentioned, these sums vanish in macroscopically large spherical systems, being nonzero otherwise. The sums 1Z {71), 
T, and V, on the other hand, contain r^ 6 (instead of r^ 3 or r ^ r Jk), which makes them rapidly convergent and 
shape-independent .0 

Figure [l] shows the thermodynamic perturbation theory susceptibility (thin lines) for — £d/2<r = 0.02 and x m 
the noninteracting limit (thick lines) for a system with parallel axes with the form of a prolate ellipsoid, a sphere, 
and an oblate ellipsoid (the symbols correspond to Monte Carlo simulations to be discussed below). For the non 
spherical systems C ^ 0,E3 and the corrections to the susceptibility are largely dominated by the first order term, 
which corresponds to the direct action discussed above. The shape dependence is easily understood by recalling 
the behavior of two dipoles: if their axes are aligned, they minimize their dipole-dipole energy lying parallel along 
the line joining them, while if this line is perpendicular to the axes, they minimize the interaction energy pointing 
along opposite directions. Therefore, in elongated systems the aligning effect dominates and \ is larger than the 
noninteracting susceptibility [Fig. |l](a)], while in oblate systems the opposite occurs with the associated decrease of 
X [Fig. 0(c)]. In the sphere the direct term exactly cancels (C = 0), and negative second order corrections (which 
incorporate indirect and reaction terms) determine the susceptibility [Fig. |l|(b)]. 

Concerning the specific heat, due to the presence of the rapidly convergent lattice sums 1Z, T, and V, this quantity 
does not depend on the shape of the system. Physically, this is a consequence of the absence of linear term in Eq. 



(3.8), which follows from (Gi) a = 0, which in turn requires the absence of a bias field. 
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B. Anisotropy dependence 



To illustrate the importance of taking the anisotropy into account, we are going to compare the thermodynamic 
perturbation theory results with: (i) those obtained by an ordinary high temperature expansion for zero anisotropy 
(the 5; — > limit of our expressions) and (ii) the results in the discrete orientation limit (Si — + 1). 

In the cases of cubic or completely disordered lattices we already know that the lattice sum V vanishes. Then, 
the first corrections to the susceptibility d ue t o the interactions become exactly independent of the anisotropy if the 
axes are distributed at random [see Eqs. (3.5)-( pV7[ ) where the only anisotropy dependent term is multiplied by V]. 
Therefore, the susceptibility coincides with that obtained by an ordinary high-temperature expansion for isotropic 
dipoles (see Ref. [t2| V = was implicitly used in that work) but also with the perturbative x f° r Ising-like spins with 
randomly distributed "Ising" axes [Fig. ^|(a)]. This generalizes to interacting spins the well known result, discussed 
in Ref. |] (see also Ref. ^Oj), i.e., the absence of anisotropy effects on the equilibrium linear susceptibility of systems 
with random anisotropy. 

Nevertheless, when V^O (as occurs in tetragonal lattices) the susceptibility will show some anisotropy dependence 
for random axes (although weak, since V ^ is accompanied by C ^ and \ is then dominated by the anisotropy 
independent first order term). In any case, on inspecting Eqs. ( |3.2| ) — ( 3.4 ) one would expect important differences for 
parallel axes between the thermodynamic perturbation theory results and those were the anisotropy is not included 
(Si — > 0). This is what actually occurs as Fig. ^(b) illustrates: the susceptibility is not only larger for parallel axes 
than for isotropic dipoles, but also the temperature dependence is stronger, owing to the extra, anisotropy-induced, 
temperature dependence of the coefficients a n via Si (a). This is clearly seen when comparing the correction terms 
Ax themselves [inset of Fig. ^(b)]. 

Note that the susceptibilities in the isotropic and Ising cases constitute lower and upper bounds to the actual X- 
The upper bound is slowly approached at low enough temperatures (a 3> 1), completing in this way the crossover from 
isotropic behavior at high T to the discrete-orientation behavior at low T. Note finally that the lowest temperatures 
displayed (a ~ 20) are still inside the range where an ordinary magnetization experiment yields the equilibrium 
response (a ~ 25). 

Concerning the specific heat, the part corresponding to the noninteracting system, & , does not depend on the 
anisotropy axes orientations. The reason is that if the spins are independent, they cannot probe the relative orienta- 
tions of their axes and there is no preferential direction to compare their orientations with (as that of the probing field 
in the susceptibility). Thus, the noninteracting specific heat (inset in Fig. [?]) only reflects the individual behavior of 
the spins in their single-spin anisotropy potentials, and its peak (at <j ~ 5) reflects the "transition" from the isotropic 
behavior at high T (with c v oc 1/T 2 ) to the discrete orientation behavior at low T. This can_be considered as a sort 
of Schottky peak, due in this case to the "depopulation" of the high-energy "barrier levels" E3 

The corrections due to the coupling depend naturally on the orientations of the axes. For a given axes distribution 
the specific heat increases with the interaction strength (Fig. ^), since an amount of heat injected in the system can 
partially be stored in the form of the potential energy of interaction. For the same reason, since the average of the 
interaction energy is larger (in magnitude) in systems with aligned axes, one would expect a larger c v in this case. 
Figure |3| shows that the effect of the interaction is indeed stronger in a system with aligned axes than in a system 
with random anisotropy. 



C. Validity limits of the perturbational results 

Since the analytic expressions derived for x and c v are expansions valid in principle for ^ < 1, which corresponds 



to T oc I/a ^> 2/i<j [Eq. (2.7)], the results will deviate appreciably from the exact quantities at sufficiently low T. In 
order to estimate the limits of validity, it would be desirable to compare the formulas with the results of some method 
treating the interactions without approximations. Besides, as the fundamentals of the expansion are the same for the 
different quantities, it would suffice to estimate the range of validity for one of them. 

We have compared the analytical x with the susceptibility computed by a Monte Carlo method (described in 
Appendix [d]), which except for statistical and finite sampling errors is exact to compute equilibrium properties. 
Returning to Fig. |l]we observe that x obtained by thermodynamic perturbation theory (thin lines) describes accurately 
the simulated susceptibility (symbols) at high and intermediate T, while the results start to deviate slightly at the 
lowest temperatures displayed (a ~ 4). Therefore, since h& = 0.02 was used in this graph, an estimate of the lower 
temperatures attainable is £d ~ 1/6, which is milder than the a priori restriction £d -C 1. 
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V. SUMMARY AND CONCLUSIONS 



We have obtained approximate analytical expressions for the equilibrium linear susceptibility and specific heat 
of classical spins interacting via dipole-dipole interactions by means of thermodynamic perturbation theory. The 
formulas account for the interactions to second order in the coupling constant over temperature but are exact in the 
anisotropy. As the results are valid for any strength and sign of the anisotropy, provided only the interaction strength 
is weak enough (£<j < 1 /6) , they include as particular cases the linear response and specific heat of isotropic as well 
as strongly anisotropic spins (discrete-orientation model and plane rotators) . 

The expressions derived also account for the different orientational distributions of anisotropy axes. For randomly 
distributed axes and sufficiently isotropic lattices (e.g., cubic or completely disordered), the linear susceptibility 
becomes independent of the anisotropy, at least to second order in the coupling constant. This extends to interacting 
systems the well known absence of anisotropy effects on the equilibrium linear response of systems with random 
anisotropy. The same holds for the corrections to the specific heat due to the interactions (indeed, without restrictions 
on the lattice type). For a general axes distribution, however, the anisotropy effects do not disappear. The importance 
of including them has been illustrated in the case of coherent axes by showing the failure of ordinary high-temperature 
expansions (for either isotropic or strongly anisotropic spins) to describe the exact susceptibility in cases where the 
thermodynamic perturbation theory yields accurate results. 
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APPENDIX A: AVERAGES WEIGHTED WITH THE UNIAXIAL ANISOTROPY BOLTZMANN 

FACTOR 

The averages we need to calculate are all of products of the form I m = (Yl™ = i(cn • s)) a , where the c n are arbitrary 
constant vectors. Introducing the polar and azimuthal angles of the spin, ($, <p), we can write I m as 

= ip 2 ^ dip Jq d$ sin ■& n» =1 (c n ■ s) exp[g(s ■ ft) 2 } 
Jq™ dip J™ dd sin i? exp [a(s-ri) 2 ] 

For odd m, J m is an integral of an odd function over a symmetric interval and hence I m = 0. To calculate the 
susceptibility and specific heat to second order in we require 1% and I4, which will be calculated using symmetry 
arguments similar to those employed to derive the a = unweighted averages (see, for instance, Ref. |2^ ). 

Note that I2 is a scalar bilinear in ci and cV The most general scalar with this property that can be constructed 
with the vectors of the problem [c\, 02, and n) has the form 

h = A ci • c 2 + B (ci • ft) (c 2 ■ n) . 

To find the coefficients A and B one chooses particular values for the c„: (i) If ci || c 2 _L n then 7 2 = A. Thus, setting 
n = z and c\ = C2 — x, one has s ■ n — cos§ = z and (ci • s)(c2 • s) = (I — z 2 ) cos 2 ip, so the integral reads 

J Q 27r dip cos 2 f J_ x dz(l — z 2 ) exp(az 2 ) 
J Q 27r dip J , dz exp(crz 2 ) 

= |[1-<A]= 1 ^ ! 

where ^(c) = (P2(z)) a is the average of the second Legendre polynomial P%{z) — i (3 z 2 — l) over the noninteracting 
distribution, (ii) If ci || c*2 || ft, then I2 = A + B. Putting n = c.\ = c*2 = z the integral is given by 

J\dzz 2 exp(az 2 ) l + 25 2 

A + B = — j = (z ) = . 

J ^ dz cxp(crz 2 ) 1 3 

Therefore, since I2 = ((ci • s){c% ■ s)) a , we get for the second order moment 







1 - s 2 

((ci • s)(c 2 ■ s)) a = — g — c i • c 2 + 5 2 (ci • n)(c 2 • ft) . 
We can similarly calculate I4 by constructing the most general scalar fulfilling certain properties, getting 

((ci • s)(c 2 • s)(c 3 • s)(c 4 • s)) a 
= A 4 [(ci • c 2 )(c 3 • c 4 ) + (ci • c 3 )(c 2 • c 4 ) + (ci • c 4 )(c 2 • c 3 )] 

+ A 2 [(ci • c 2 )(c 3 • ra)(c 4 • ra) + (ci • c 3 )(c 2 • n)(c 4 • ra) 

+ (ci • c 4 )(c 2 • r?)(c 3 • ft) + (c 2 • c 3 )(ci • ft)(c 4 ■ ft) 
+ (c 2 • c 4 )(ci • n)(c 3 • n) + (c 3 • c 4 )(ci • n)(c 2 • n)] 

+ 5 4 (c*i • n)(c 2 • n)(c 3 • n)(c 4 • n), 

where A 2 and A 4 are combinations of the first Si (a) 



(Al) 



A 2 = i(5 2 -5 4 ), 



A 4 = ^- 



25 2 



35 21 



1 

15 



(A2) 



(A3) 



Therefore, Eq. (A2) involves S2 as well as Si(a) = (P 4 (z)) a , the average of the fourth Legendre polynomial Pi(z) 
I (35 z A - 30 z 2 + 3) with respect to W a . 

Finally, introducing the following tensor and scalar shorthands 



1-52 
r = — - — 1 + S 2 ft n , 

nr- A 2 A 2 
A = VAll + -^=nn , = 5 4 -3-^, 
VA 4 A 4 

where 1 is the identity tensor, the results for the moments can compactly be written as 

((ci ■ s)(c 2 ■ s)) a =(c 1 -r-c 2 ) 
((ci • s)(c 2 • s)(c 3 • s)(c 4 • s)) a = (ci • A • c 2 )(c 3 • A • c 4 ) 

+ (ci • A • c 3 )(c 2 • A • c 4 ) 
+ (ci • A - c 4 )(c 2 • A • c 3 ) 
+ fi(ci • n)(c 2 • n)(c 3 • n)(c 4 ■ ra) , 

which facilitates the manipulation of the observables. 

The quantities Si can be computed using the following homogeneous three-term recurrence relatior 



1 



•la 



(21 - 1)(2/ + 3) 



2a 
21 + 1 



I - 1 
21-1 



Si- 



knowing the first two terms: Sq = 1 and S 2 , which is given by 



1 

2^ 



1-2 



1 + 2 
21 + 3 



S, 



1+2 



, 



(A4) 
(A5) 

(A6) 

(A7) 

(A8) 
(A9) 



The one-spin partition function Z a = J dz exp(crz 2 ) can be written in terms of error functions of real and "imagi- 
nary" argument as 



'n/a erfi(- s /o r ) , a > 



VVHerf(vVf), <J<0 



(A10) 



The less familiar erfi(a;) is related with the Dawson integral D(x), so in the easy-axis case one can write Z a = 
(2e CT / < y r a)D( y /a) and compute D(x) with the subroutine DAWSON of Ref. ||. 

Note finafe-that in the isotropic limit (Si — > 0), Eqs. (Al) and (A2) reduce to the known moments for the isotropic 
distributior£§B 



((ci • s)(c 2 ■ s)) iso = i ci • c 2 , 
((ci ■ s)(c 2 ■ s)(c 3 - s)(c 4 ■ s))iso = jg[(ci ■ c 2 )(c 3 • c 4 ) 

+ (ci • c 3 )(c 2 • c 4 ) 
+ (ci • c 4 )(c 2 • c 3 )] 



(All) 



(A12) 
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These expressions are formally identical to those for the average of a quantity involving the anisotropy axes n,-, when 
these are distributed at random Y2i f(ni) — > J = /■ For instance, for arbitrary n-independent vectors v\ 

and v 2 , we have 

■ft ■ ni)(v 2 ■ Hi) — * (vi ■ n)(v 2 ■ ft) = § v x ■ v 2 . (A13) 



APPENDIX B: GENERAL FORMULAE FOR THE COEFFICIENTS OF THE SUSCEPTIBILITY 



The general expression for the equilibrium linear susceptibility is given by Eq. (3.1) with the following expressions 
for the coefficients 



a = 



a i 



"2 



— V h ■ Ti ■ h 
N ^ 

I 



N 



i jjti 



i jjti k^j 



■{yy. 1 



(h- Ai ■ h)(f i3 ■ A; ■ fi 



+ 2(h-Ai-hjY 



+ S 2 



+ Q(h ■ fii) 2 (ni ■ hjY 
(h- Ai - h)(nj ■ Gij ■ Ai ■ Gy ■ Hj) 
+ 2(h ■ Ai ■ ■ Hj) 2 
+ Q(h ■ ni) 2 (fii ■ Gij ■ fij) 2 
l-S 2 



i jjii 



-{fij ■ Ti ■ fij) 



S 2 (nj ■ Gij ■ Ti ■ Gj 



where Gij, fij and fij arc defined in Eq. (|2.3| ), and T, A, and fl in Eqs. (|A4| )- (|A5[) and also involve the Si(a). 

When calculating these coefficients, the same type of averages appear as in the isotropic case (see Refs. mJ13 for 
details of the calculation) and with the same multiplicities. The only difference is the weight function and hence the 
formulas required to calculate those averages [Eqs. (fO) and (pO) instead of Eqs. ([All]) and flA12])]. 



APPENDIX C: GENERAL FORMULA FOR THE COEFFICIENT B 2 OF THE SPECIFIC HEAT 



In the general expression fl3.9| ) for the specific heat the coefficient bo is given by Eq. (3.1C), while b 2 reads 
Nb 2 = |{2(1 - S 2 ) 4aS 2 a 2 S' 2 '} ]T £ 

i j^i 

+i{2S 2 (l - S 2 ) + iaS' 2 (l - 2S 2 ) 
+ a 2 [S' 2 l {l-2S 2 )-2{S' 2 ) 2 ]} 
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i j^i 

+ {S 2 2 + AaS 2 S' 2 + a 2 [S 2 S' 2 ' + (S' 2 ) 2 }} J2J2^ ' G « ' H ^ > ( C1 ) 

where /' = df /da. An arbitrary 5; can be expressed in terms of averaged Legendre polynomial by means of the 
following differential-recurrence relation 

dSi _ (I- 1)1 21(1 + 1) 

da (2Z — 1)(2Z + 1) l ~ 2 3(2Z - 1)(2Z + 3) ' 



(2l + l)(2l + 3) ^ 3 

This useful formula can readily be demonstrated by taking the derivative of the definition of Si = (Pi 

f\dzPie az2 / f\dze az2 : 

dSi J\dz z 2 P l e°* 2 J^dz Pi e az2 j\dz z 2 e 



d<? f\dze™ 2 (j\dze^) 2 

= (z 2 P) a - |S, (1 + 2S 2 ) , (C3) 

where (z 2 ) & = (1 + 2S 2 )/^ has been used. The product of Pi with z can be expanded in Legendre polynomials by 
using the corresponding relation for associated Legendre functions P™ (Ref. Ch. 12) 

*pi = [ipi-i + + m+i] • (C4) 



Multiplying zPi by z, using again Eq. (C4) to expand zPi±\ on the right-hand side, and gathering terms yields 



Z 2 R _ (l-W p , 2l(l + l)-l 
1 (2/ — 1)(2Z + 1) (2/ - 1)(2Z + 3) 

(/ + !)(/ + 2) 



-P 



1+2 



(2Z + 1)(2Z + 3) 

Averaging this result and substituting in Eq. (jC3|), we finally get the desired Eq. (p2| 



APPENDIX D: MONTE CARLO SIMULATIONS 

In order to obtain nonperturbative results to test the analytical expressions we have performed careful Monte Carlo 
simulations, by a method similar to that employed in Ref. The trial Monte Carlo step (MCS) is a random rotation 
of the spin within a cone, achieved by generating a random vector with uniform probability on the surface of a sphere 
of radius g, adding the random vector to the initial spin, and finally normalizing the resulting vector. For a better 
acceptance ratio of the generated configurations, we scale that radius with T as g = 0.7 /y/a. 

In terms of the maximum anisotropy field hqHk = 2A/m, where A is the anisotropy parameter and m is the 
magnetic moment, a dimensionless probing field can be defined as AH/Hk or, in temperature units, as A£ = 
2aAH/H K . The latter is the argument in the Langevin function for the (isotropic) magnetization, and it controls 
if the response is linear in AH. Therefore, to treat all the temperatures on the same footing A£ oc AH/T is kept 
constant in the simulation (A£ = 0.15) which requires to decrease the probing field with T. The value A£ = 0.15 
used should ensure linear response since it is well below the values A£ = 0.3-0.5, where nonlinear terms can start to 
contribute appreciably to the equilibrium response (see, for instance, Appendix C of Ref. p7| ). 

The number of dipoles in each simulation is indicated in the panels of Fig. u\ an d the simulations are done with 
open bo unda ry conditions. The absolute value of the coefficient in front of Eq. (2.2) can be written as 2 A h^, so that 
hd [Eq. ( |2.7| )1 is used as the input in the simulation. In order to control that the response is the equilibrium one, we 
apply a sinusoidal probing field of low frequency (/ = 8 x 10~ 6 MCS -1 ) and check that further reducing / does not 
change the results. Besides, the 10 first periods are excluded from the average and the susceptibility sampled over 
the 20 following periods. 
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Prolate spheroid: a=b=2, c=6; N=97 




Oblate spheroid: a=b=4, c=2; N=125 



FIG. 1. 
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Equilibrium linear susceptibility 




0.35 0.4 0.45 0.5 
1/0 oc T 

temperature for three different ellipsoidal systems with equation 



x /a + y /b + z /c < 1 resulting in a system of N dipoles. The susceptibility is given in reduced units \ — xiHic/m), the 
spatial arrangement of the spins is simple cubic and the probing field is applied along the anisotropy axes, which are parallel to 
the 2-axis. The thick lines are the equilibrium susceptibility of the corresponding noninteracting systems (they are equal in all 
three cases); thin lines are the s usce ptibilities including the corrections due to the dipolar interactions obtained by thermody- 
namic perturbation theory [Eq. (3.1)]; and the symbols represent the susceptibility obtained with a Monte Carlo method. The 
dipolar interaction strength is hd = £d/2<7 = 0.02. (For the prolate and oblate ellipsoids the second order correction is very 
small in the temperature interval displayed and omitting it the curves visually coincide.) 
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Random anisotropy axes and cubic lattice. Prolate spheriod 




Parallel anisotropy axes and cubic lattice. Prolate spheriod. 
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FIG. 2. Equilibrium linear susceptibility vs. temperature for the same prolate ellipsoid as in Fig. [l](a) and the dipolar 
interaction strength hd = £d/2<7 = 0.004. The spins are arranged on a simple cubic lattice with: a) randomly distributed 
anisotropy axes and b) parallel anisotropy axes. Thick lines are the susceptibilities for independent spins and thin lines are the 
susceptibilities obtained by thermodynamic perturbation theory. For comparison we have displayed x obtained by a classical 
high temperature expansion for isotropic spins (crosses) and for Ising spins (dashed line). Inset: Comparison of the corrections 
to the susceptibility due to interactions (Ax = X~ Xnon-int). (Note that the temperatures displayed are above the lower validity 
limit £d ~ 1/6 estimated in the text: 1/a — 2ftd/£d ~ 0.048.) 
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FIG. 3. Specific heat per spin vs. temperature for non-interacting spins (thick line), and weakly interacting spins with 
randomly distributed anisotropy axes (dashed lines) and parallel axes (thin lines) arranged on a simple cubic lattice. In each 
case, hd = £d/2cr = 0.003 and 0.006 from bottom to top. The inset shows the specific heat for non- interacting spins over a 
wider temperature interval. 
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